A multi-omics approach to elucidate the mechanisms of action of a dietary muramidase administered to broiler chickens

A novel dietary muramidase has been shown to have positive effects on broiler chickens. However, very little is known about its mechanisms of action. The present multi-omics investigation sought to address this knowledge gap. A total of 2,340 day-old male broilers were assigned to 3 groups (12 replicates each) fed, from 0 to 42 d, a basal diet (control group—CON) or the basal diet supplemented with muramidase at 25,000 (low-dose group—MUL) or 45,000 LSU(F)/kg feed (high-dose group—MUH). MUH significantly outperformed CON in terms of cumulative feed intake (4,798 vs 4,705 g), body weight (2,906 vs 2,775 g), and feed conversion ratio (1.686 vs 1.729), while MUL exhibited intermediate performance. At caecal level, MUH showed the lowest alpha diversity, a significantly different beta diversity, a reduction in Firmicutes, and a rise in Bacteroidetes, especially compared with MUL. MUH also exhibited a considerable decrease in Clostridiaceae and an overrepresentation of Bacteroidaceae and Lactobacillaceae. At blood level, MUH had lower hypoxanthine—probably due to its drop at caecal level—histidine, and uracil, while greater pyruvate, 2-oxoglutarate, and glucose. This study sheds light on the mode of action of this muramidase and lays the groundwork for future investigations on its effects on the intestinal ecosystem and systemic metabolism of broiler chickens.


Results
Performance traits. Chicks weighed approximately 42 g at placement with no inter-group significant differences. At the end of starter phase, MUL and MUH exhibited a higher BW than CON (199.5,204.8, and 205.8 g for CON, MUL, and MUH, respectively; p < 0.05), whereas only MUH showed a lower FCR than CON (1.267 vs 1.240 for CON and MUH, respectively; p < 0.05) ( Table 1). MUH reached the greatest BW at the conclusion of the first grower phase (765.3, 785.1, and 818.8 g for CON, MUL, and MUH, respectively; p < 0.05) as well as a higher FI than CON (866.0 and 907.7 g for CON and MUH, respectively; p < 0.05) ( Table 1). At the end of the second grower phase, MUH showed the greatest BW (1,344, 1,375, and 1,443 g for CON, MUL, and MUH, respectively; p < 0.05) and lowest FCR (1.673, 1.651, and 1.590 for CON, MUL, and MUH, respectively; p < 0.05) ( Table 1). At the conclusion of finisher phase, MUH exhibited a greater BW than CON (2,775 and 2,906 g for CON and MUH, respectively; p < 0.05) and MUL (2,835 and 2,906 g for MUL and MUH, respectively; p = 0.05), while other performance traits were unaffected ( Table 1). Results of the overall trial indicate that, in addition to BW, MUH outperformed CON in terms of cumulative FI and FCR (4,705 g and 1.729, and 4,798 g and 1.686 for CON and MUH, respectively; p < 0.05) ( Table 1; Fig. 1). Polynomial contrasts also revealed that cumulative FI and BW significantly increased while FCR decreased across groups in a linear fashion (Fig. 1). Lastly, mortality rate was not significantly influenced (Table 1).
Processing yields. Data of carcass and cut-up yields-generated at processing in a commercial plantwere not subjected to statistical analysis because measured on a group basis. Notwithstanding, the considerable sample size (i.e., more than 750 observations/group) made it possible to notice that muramidase-supplemented groups had a greater eviscerated carcass yield (70.1, 70.4, and 70.8% for CON, MUL, and MUH, respectively) and breast yield calculated as percentage of eviscerated carcass weight (30.6, 30.9, and 31.3% for CON, MUL, and MUH, respectively).
Muramic acid concentration in excreta samples. Concentration of muramic acid (total, soluble, and their ratio), used as a measure for the concentration of hydrolyzed bacterial PGN in excreta samples, is given in Table 2. In every feeding phase, muramidase-supplemented groups exhibited a significantly higher soluble fraction of muramic acid as well as a greater ratio with respect to total muramic acid than CON. Moreover, it was found a weak negative Pearson's correlation coefficient between soluble to total muramic acid ratio and FCR (r = −0.30; p < 0.001).
Incidence and severity of foot-pad dermatitis and breast muscle myopathies. Occurrence of foot-pad dermatitis (FPD) was significantly associated with the factor group (Fig. 2). MUH was 0.58 times less likely to develop FPD than MUL; that is, supplementing chickens with muramidase at high dose decreased by 42% the relative risk of developing FPD compared to their counterparts supplemented at low dose. However, MUL diet tended to increase by 33% the relative risk of FPD development compared to CON diet (Table 3). On the other hand, breast muscle myopathies did not show a significant relationship with the factor group (Fig. 3). www.nature.com/scientificreports/ Caecal microbiome. MUH exhibited a lower alpha diversity than MUL (p < 0.05 for Simpson and p = 0.06 for Shannon and Inverse Simpson, respectively), while CON did not differ from muramidase-supplemented groups (Fig. 4). Regarding beta diversity, the PCoA revealed an evident segregation of MUH samples (Fig. 5).
Besides this visual distinction, the PERMANOVA confirmed a group effect on beta diversity (p = 0.005, R 2 = 0. 22 (Fig. 7). Lastly, at species level, relative abundance of C. phytofermentans, C. saccharolyticum, C. cellulolyticum, and C. butyricum was lower in MUH than other groups (Fig. 8). This variation also applies to Eubacterium rectale, Roseburia intestinalis, Ruminococcus albus, C. perfringens, C. botulinum, and Listeria monocytogenes (Fig. 8-9). An www.nature.com/scientificreports/ opposite pattern, however, was observed for Bacteroides thetaiotaomicron whose relative abundance was higher in MUH than MUL (p < 0.05) (Fig. 8). Relative abundance of genes that significantly differed between groups is shown in Fig. 10. At glycan biosynthesis and metabolism level, MUH showed a greater relative abundance of genes associated to glycosaminoglycan degradation pathway than MUL (p < 0.02), while relative abundance of genes involved in peptidoglycan biosynthesis pathway was affected in the opposite fashion (p < 0.02). Genes involved in starch and sucrose metabolism and amino sugar and nucleotide metabolism were affected by the factor group: the former had a higher relative abundance in CON and MUL than MUH (p < 0.01), whereas the latter showed an increasing trend in MUL compared to MUH (p < 0.1). Genes involved in seleno compound metabolism pathway showed an increase in MUL compared to MUH (p < 0.05). However, MUH exhibited a higher   www.nature.com/scientificreports/ relative abundance of genes involved in glutathione metabolism pathway than both CON and MUL (p < 0.02 and p < 0.05, respectively). Lastly, two transport and catabolism pathways revealed an inter-group difference: genes involved in lysosome path had a greater relative abundance in MUH than MUL (p < 0.05), while genes linked to peroxisome path varied the other way around (p < 0.05).
Plasma and caecal content metabolomes. Plasma and caecal content 1 H-NMR spectra were registered and 54 and 78 molecules were assigned and quantified, respectively. At caecal content level, the concentration of 4 metabolites showed a significant inter-group difference (Table 4). While acetate, ferulate, and formate were greater in MUL than MUH (p < 0.05), hypoxanthine was higher in CON than MUH (p < 0.05). The rPCA model shown in Fig. 11 was built on these molecules. The principal component one (PC1) accounts for 68% of the variance explained by the model and summarizes the differences between groups. PC1 scores of MUH samples are markedly or marginally higher than those of other groups, resulting in a group-based clustering of samples mainly led by ferulate and formate (r < −0.5). At plasma level, the concentration of 9 metabolites showed a significant variation between groups (Table 5). Specifically, a higher concentration of pyruvate was observed in MUH than CON and MUL (p < 0.05), while 2-oxoglutarate, glucose, and uridine were greater in MUH than MUL (p < 0.05). On the other hand, MUH showed a higher concentration of myo-inositol than CON (p < 0.05), whereas CON exhibited a greater concentration of histidine than MUH (p < 0.05). Similarly, both CON and MUL had a higher concentration of hypoxanthine than MUH (p < 0.05), while MUL showed a greater uracil concentration than MUH (p < 0.05). Figure 12 illustrates the rPCA model produced as described above. Samples of MUH are characterized by higher PC1 scores and distinctly segregate (p < 0.05): this separation is predominantly driven by pyruvate, 2-oxoglutarate, glucose, uracil, and hypoxanthine (r > 0.5 or < −0.5).

Figure 2.
FPD incidence and severity (score 0, no lesions; score 1, mild lesions; score 2, severe lesions) of CON, MUL, and MUH at slaughter (42 d). FPD were macroscopically measured on 1 foot per bird (n = 675, 673, and 684 for CON, MUL, and MUH, respectively). Table 3. FPD risk ratio computation on 2 by 2 tables aligning the combinations of levels of the factor group and putting binarily aggregated FPD scores in columns. † The inverse contrasts produce a risk ratio equals to the reciprocal of the risk ratio shown. ‡ 95% Wald confidence interval is given in brackets. § Calculated as risk ratio minus 1 percentagewise. ¤ Count data were statistically analyzed via Pearson's chi-square test. www.nature.com/scientificreports/

Discussion
The purpose of this multi-omics investigation was to shed light on the mechanisms of action of a dietary muramidase supplemented to broiler chickens. MUH, the experimental group receiving the muramidase at high dose (i.e., 45,000 LSU(F)/kg feed), significantly outperformed the control group CON in terms of cumulative FI, BW, and FCR at 42 d. On the other hand, MUL, the low-dose group (i.e., 25,000 LSU(F)/kg feed), showed intermediate cumulative performance and did not differ from CON in a significant manner. It is worth highlighting that cumulative FI, BW, and FCR improved proportionately with muramidase dose. These results broadly support those of previous research assessing the administration of the same muramidase to broiler chickens, wherein birds supplemented at high inclusion levels (i.e., 35-45,000 LSU(F)/kg feed) performed better than their control and low-dose counterparts [21][22][23][24] . The non-invasive technique employed here to measure the muramidase-mediated PGN hydrolysis in excreta samples can be a reliable alternative to the ex vivo analysis lately illustrated by  www.nature.com/scientificreports/ Frederiksen and co-workers 25 . Indeed, the proposed method confirmed that this muramidase effectively hydrolyzes bacterial PGN causing the release of fragments that, according to recent reports, can support intestinal health and performance of broiler chickens [22][23][24] . The observed reduction in FPD occurrence associated with muramidase supplementation, especially at high dose, is consistent with results obtained by Pirgozliev and colleagues 24 . Poultry foot-pad welfare greatly depends on litter quality (e.g., moisture and ammonia levels) and its management 26 . Although measuring litter parameters was beyond the scope of this work, a better nutrient utilization as well as less watery excreta-both commonly resulting from enhanced feed efficiency-may have played key roles in FDP risk attenuation 27 . Analysis of breast fillets revealed that the muramidase did not affect WS, WB, and SM. Therefore, under our experimental settings, this muramidase was able to improve growth performance without exerting negative effects on the occurrence of breast muscle myopathies, with positive implications for the sustainability of poultry meat production. Feeding trials on pigs 10,28-31 , rabbits 12 , and broiler chickens [13][14][15] suggest that the beneficial effects of dietary muramidases can be ascribed to a modulatory activity  www.nature.com/scientificreports/ on the GI microbiota. In the present investigation, the muramidase supplemented at high dose produced a drop in caecal alpha diversity. This result is in agreement with those of previous studies on pigs 29 and, above all, on the same supplement fed to broiler chickens 19,23 . Moreover, MUH exhibited a different bacterial community structure at caecal level, especially compared to MUL. This is in accord with research on piglets and lactating sows 29,30 , and supports findings obtained in broiler chickens treated with the same muramidase 19 . Not only was changed the overall caecal bacterial community structure, but also its taxonomic composition. The observed underrepresentation of Firmicutes and outgrowth of Bacteroidetes, particularly evident for the comparison between MUH and MUL, seem to be consistent with results of Maga and colleagues 28 . These researchers fed pigs with the milk produced by transgenic goats expressing the human muramidase and proved that, in fecal samples of treated pigs, abundance of Firmicutes fell whereas that of Bacteroidetes raised over time. In the current study, MUH also showed a significant decrease in the Firmicutes to Bacteroidetes ratio compared with MUL. A considerable amount of papers has been published on the role played by the Firmicutes to Bacteroidetes ratio in the microbiota-to-host energy supply and development of obesity. However, since contradictory outcomes are not uncommon, Magne et al. 32 have advocated that a direct causality between this ratio and health status of the host is hard to be attested. A possible explanation for the detected differences in abundance of Firmicutes and Bacteroidetes is that bacteria positive to the Gram staining, like Firmicutes, are generally more vulnerable to the hydrolytic action of muramidases on PGN. Indeed, Firmicutes possess an undefended, thicker cell wall lacking  www.nature.com/scientificreports/ an outer protecting lipid membrane and offering up to 40 PGN layers as a substrate to muramidases 8 . The lower Firmicutes abundance showed by MUH may explain the observed decrease in genes associated with peptidoglycan biosynthesis pathway. Moreover, the reduction in genes involved in amino sugar and nucleotide metabolism pathway, in which NAG is directly involved 33 , can be taken as another indicator for the inhibition of bacteria with high PGN synthesis capacity. MUH showed a considerable drop in Clostridiaceae. Interestingly, earlier studies on pigs 28 and rabbits 12 established that dietary supplementation of muramidases causes a depression in GI Clostridia. Furthermore, Sais et al. 23 found a decreasing trend in Clostridium count at ileal level after supplementing broiler chickens with the muramidase tested here, yet at 35,000 LSU(F)/kg feed. MUH also exhibited a decrease in several butyrate-producing Clostridia, Lachnospiraceae (viz., Roseburia intestinalis and Ruminococcus albus), and Eubacterium rectale 34 . This was contrary to expectations as butyrogenic bacteria have traditionally been linked to gut health 35-37 , while chickens have been shown to benefit from short-chain fatty acids, especially butyrate, released by GI bacteria [38][39][40][41] . Future studies on this topic are therefore recommended. C. perfringens is the causative agent of necrotic enteritis, a gut disorder that causes to the poultry sector a financial burden of 2 billion dollars yearly 42 . Hence, the control of C. perfringens is vital, especially in the antibiotic growth promoters-and antimicrobial-free era. Even though the abundances found in this study were low, MUH showed a reduction in C. perfringens. This finding is consistent with that of Liu et al. 13 who hindered the intestinal colonization of C. perfringens in broiler chickens orally challenged with this pathogen and treated with a dietary muramidase. Poultry are also susceptible to C. botulinum neurotoxins and can sporadically manifest avian botulism, a flaccid paralytic disease 43 . Despite the low abundances detected here, the decrease in C. botulinum exhibited by MUH is an issue that deserves further research. L. monocytogenes, an important human pathogen 44 , showed a minor presence in MUH although the measured abundances were rather low. Field studies revealed that poultry can be a reservoir of L. monocytogenes, thereby contributing to the contamination of processing facilities 45 . Interestingly, the observed inhibition of C. botulinum and L. monocytogenes supports earlier studies indicating that muramidases are effective solutions against these pathogenicbacteria 46,47 . The higher abundance of Bacteroidaceae in MUH agrees with results obtained in pigs 28,30 . However, contrary findings have also been found in the latter species 29 . Bacteroides have positively been associated to human gut health due to their propionate-producing ability 48,49 , while B. thetaiotaomicron has been included in a probiotic blend to restore gut eubiosis after antibiotic therapies 50 . Therefore, it can conceivably be posited that caecal Bacteroidaceae can promote intestinal health of broiler chickens as well. The observed increase in Lactobacillaceae is comparable to previous results on both non-avian species, like pigs 28-30 and rabbits 12 , and broiler chickens supplemented with the muramidase used here 21,23 . The rise in Lactobacillaceae also differs from previous findings 10, 19 and can be  www.nature.com/scientificreports/  Figure 11. rPCA model on caecal metabolites of Table 4. In the score plot (a), samples of CON ("A"), MUL ("B"), and MUH ("C") are drawn as squares, circles, and triangles, respectively. Wide circles are the group medians. The box plot (b) summarizes the position of samples along PC1. The loading plot (c) reports the correlations between the concentration of each metabolite and its importance over PC1. Grey bars indicate significant correlations (p-value < 0.05). www.nature.com/scientificreports/ contradictory considering the abovementioned proneness of Gram-positive bacteria to muramidase-mediated PGN hydrolysis. However, some lactic acid bacteria employed in the production of hard-cheese have been shown to be or gradually become resistant to muramidases 51,52 . Several Lactobacillaceae have also probiotic features 53 that might have supported the performance of muramidase-supplemented birds, especially those treated at high dose. Attributing to this muramidase stimulating or, at least, non-inhibiting effects upon enteral Lactobacillaceae warrants further investigations. Results of KEGG analysis can help interpret the changes occurred in the caecal metabolome. In MUH, the decrease in abundance of genes of starch and sucrose metabolism path can be behind the measured reduction in fermentations-deriving organic acids, such as acetate, ferulate, and formate. The enrichment in genes linked to glutathione metabolism path-a bacterial cells' antioxidant tool 54 -can be associated with hypoxanthine drop at caecal level. Hypoxanthine is a noxious end-product of purine-catabolism, considered as a biomarker for oxidative stress [55][56][57] . Therefore, the lower concentration of hypoxanthine may have positively influenced the GI ecosystem of MUH birds. In addition, hypoxanthine decrease at intestinal level might have been the reason for its minor presence at plasma level. Lower circulating hypoxanthine, and histidine and uracil-protein-and nucleotide-catabolism end-products, respectively 56,57 -indicate that the degradation of proteins and nucleotides, intended to generate energy, may have occurred to a lesser extent in MUH birds. The higher abundance of caecal Bacteroidaceae may have increased the supply of propionate for the hepatic gluconeogenesis in MUH birds 58 , thereby leading to the observed rise in plasmatic energetic compounds such as pyruvate, 2-oxoglutarate, and glucose. The enrichment in bioenergetic compounds and reduction in prooxidant protein-and nucleotidecatabolites suggest that a more balanced energy metabolism may have stimulated the performance of high-dose supplemented birds. Surprisingly, cumulative FI, BW, and FCR were influenced in a dose-dependent fashion, while muramidase-supplemented groups showed the most marked microbiome and metabolome divergences. A possible explanation for this is that the cumulative performance benefited from an additive effect of each feeding phase, whereas the molecular outcomes at slaughter cannot fully justify the GI and metabolic dynamics of the entire grow-out period. Taken together, these findings contribute in several ways to our understanding of the mode of action of this dietary muramidase. The present study also lays the groundwork for future investigations on the effects of this muramidase on the GI ecosystem and systemic metabolism of broiler chickens. www.nature.com/scientificreports/

Methods
Experimental design, housing, and husbandry conditions. A total of 2,340 day-old male Ross 308 broilers, obtained from the same breeder flock and hatching batch, were provided by a commercial hatchery. After hatch, they were vaccinated against infectious bronchitis virus, Marek's disease virus, Newcastle and Gumboro diseases, and coccidiosis. Birds were housed in an experimental poultry facility, and randomly assigned to 3 groups (12 replicates/group) fed a commercial corn-wheat-soybean basal diet (control-CON) or the basal diet supplemented with a dietary muramidase (Balancius®, DSM Nutritional Products) at 25,000 (low-dose group-MUL) or 45,000 LSU(F)/kg feed (high-dose group-MUH) for the entire trial (0-42 d). Table 6 reports the basal diet formulation according to the 4-phase feeding program used (0-9 d, starter; 10-21 d, grower I; 22-28 d, grower II; 29-42 d, finisher). For each feeding phase, the mash basal diet was part of the same batch, while the powdery additive was added on top. The analytical inclusion levels of muramidase met the abovementioned targets. Each replicate was assigned to one of 36 floor pens (5.9 m 2 /pen) arranged in a randomized complete block design. Pens were equipped with two feeders, nipple drinkers, and chopped straw as bedding. Birds were manually fed and watered ad libitum on a daily basis. At each feeding phase switch-uniformly performed for all groups-feeders were emptied, cleaned, and refilled, while residuals weighed. The environmental temperature was modified according to the flock age by following the breeding company's recommendations. The artificial photoperiod was 23L:1D during the first 7 and last 3 d, while 18L:6D for the remainder days. Birds were han- Table 6. Basal diet composition according to feeding phases. † The premix included at 0.1% provides 1,000 FTU per kg of feed. ‡ The premix provides the following per kg of feed: vitamin A (retinyl acetate), 13,000 IU; vitamin D3 (cholecalciferol), 4,000 IU; vitamin E (DL-α_tocopheryl acetate), 80 IU; vitamin K (menadione sodium bisulfite), 3 mg; riboflavin, 6.0 mg; pantothenic acid, 6.0 mg; niacin, 20 mg; pyridoxine, 2 mg; folic acid, 0.5 mg; biotin, 0.10 mg; thiamine, 2.5 mg; vitamin B 12  Performance traits measurement. On a replicate basis, the number and body weight (BW) of birds were recorded at housing (0 d), each feeding phase switch (10, 22, and 29 d), and slaughter (42 d), while feed intake (FI) for each feeding phase. Daily weight gain (DWG), daily feed intake (DFI), and feed conversion ratio (FCR) were calculated for each feeding phase and the whole rearing period (0-42 d). The number and BW of dead or culled birds were considered to correct performance data for mortality.
Excreta collection and PGN hydrolysis assay. At the end of each feeding phase, fresh excreta samples were collected on a replicate basis (12 specimens/group; 36 specimens/feeding phase) and analyzed to evaluate bacterial PGN hydrolysis. The freeze-dried samples were resuspended and centrifuged. While the supernatant contained soluble PGN, the precipitate was enriched in insoluble PGN. Later, samples were subjected to acid hydrolysis to measure total and soluble PGN through liquid chromatography-mass spectrometric quantification of muramic acid used as a marker for hydrolyzed PGN (Novozymes A/S Biologiens, Lyngby, Denmark). Insoluble fraction of PGN was calculated as the difference between total and soluble PGN amount expressed as muramic acid. The tested muramidase has been shown to hydrolyze PGN of bacterial debris both in vitro and in ex vivo digesta samples of broiler chickens 25 . Therefore, the above described assay was performed to test, via a non-invasive method, the hypothesis of a larger proportion of hydrolyzed PGN in excreta of muramidasetreated birds.
Processing traits, breast muscle myopathies, and foot-pad dermatitis evaluation. At slaughter (42 d) in a commercial plant, groups were clearly identified and separately processed. On a group basis, carcass and cut-up yields were measured on all processed birds according to standard commercial procedures.
Plasma and caecal content collection. From

Plasma and caecal content 1 H-NMR analysis. An 1 H-NMR analysis solution with D 2 O, containing
3-(trimethylsilyl)-propionic-2,2,3,3-d4 acid sodium salt (TSP) 10 mmol/L and NaN 3 2 mmol/L was made. The solution pH was set at 7.00 ± 0.02 by phosphate buffer 1 M. TSP was selected as a reference for NMR chemicalshift, while NaN 3 was employed to avoid microorganism proliferation. Plasma samples were prepared for 1 H-NMR analysis by centrifuging 1 mL of each sample (18,630 g; 900 s; 4 °C). 0.7 mL of supernatant was subsequently added to 0.1 mL of the 1 H-NMR solution. Lastly, each sample was centrifuged once again as described above. Likewise, caecal content was prepared by vortex mixing approximately 80 mg of each sample with 1 mL of bi-distilled water: 0.7 mL of supernatant was treated as previously described for plasma. 1 H-NMR spectra were registered (600.13 MHz; 298 K) with an AVANCE™ III spectrometer (Bruker, Milan, Italy), equipped with Topspin v3.5 software. We suppressed signals from broad resonances due to large molecules with a CPMG-filter composed by 400 echoes with a τ of 400 µs and a 180° pulse of 24 µs, for a total filter of 330 ms. The water residual signal was suppressed by means of presaturation. This was done by employing the cpmgpr1d sequence, part of the standard pulse sequence library. Each spectrum was acquired by summing up 256 transients constituted by 32,000 data points encompassing a window of 7184 Hz, separated by a relaxation delay of 5 s. www.nature.com/scientificreports/ 1 H-NMR spectra were phase-adjusted in Topspin v3.5 and then exported to ASCII format by means of the built-in script convbin2asc. Spectra were processed with R software (R Core Team, 2020) through in-house developed scripts. We baseline-adjusted spectra by distinguishing baseline imperfection from NMR signals according to the "rolling ball" principle 69 implemented in the R package baseline 70 .
Signal assignment was performed by comparing their chemical shift and multiplicity with Human Metabolome Database 71 and Chenomx software library (Chenomx Inc., Edmonton, Canada, v10).
Plasma molecule concentrations were assessed by quantifying the molecules of the first sample analyzed by means of an external standard. Differences in water content between samples were then taken into consideration by probabilistic quotient normalization (PQN) 72 . Molecule concentrations in caecal samples were assessed as described for plasma by considering as reference the sample with the median water dilution assessed via PQN. The quantification of each molecule was performed through rectangular integration, considering one of its signals free from interferences. Statistical analysis. Performance data were analyzed through a one-way ANOVA with blocks, group as the experimental factor, and replicate as the experimental unit. Tukey's honestly significant difference test was used to separate the groups' means. Furthermore, polynomial contrasts were carried out to test for linear and quadratic trends in overall performance data (0-42 d). Data of carcass and cut-up yields were not subjected to statistical analysis because measured on a group basis. For each feeding phase, differences in soluble muramic acid, total muramic acid, and the ratio between them were analyzed through one-way ANOVA as described above. Pearson's correlation coefficient between soluble to total muramic acid ratio and FCR data from each feeding phase was computed and tested for significance. Count data of FPD, WS, WB, and SM were analyzed by means of Pearson's chi-square test and Fisher's exact test involving all groups and using the sampled animal as experimental unit. All these analyses were performed in R 67 with a significance level of 0.05. Moreover, FPD count data were arranged in 2 by 2 contingency tables aligning a combination of levels of the factor group (i.e., CON and MUL; CON and MUH; MUL and MUH) and having binarily aggregated FPD scores in columns (i.e., "FPD presence" as a sum of score 1 and score 2 counts; "FPD absence" as score 0 counts). Incidence risk ratio was computed on these 2 by 2 tables with epiR 73 package of R 67 . If incidence risk ratio was significant at 95% confidence interval, the risk of developing FPD was computed as incidence risk ratio minus 1 and expressed as percentage. Ecological diversity indices were analyzed at genus level with vegan 74 package of R 67 . Shannon, Simpson, and Inverse Simpson indices were chosen for alpha diversity, while the Bray-Curtis distance matrix method for beta diversity analysis. Alpha indices were analyzed with a one-way ANOVA and Tukey's post-hoc test by considering the group as experimental factor and each sampled animal as experimental unit. Beta diversity was graphically explored through principal coordinates analysis (PCoA), and analyzed with PERMANOVA-"adonis" procedure with 10,000 permutations-followed by pairwise permutation MANOVA with RVAideMemoire 75 package of R 67 . The matrix of caecal bacteria abundances was normalized for total read number in each sample and analyzed in STAMP v2.1.3 68 by using Kruskal-Wallis H-test and Games-Howell post-hoc test with group as experimental factor and each sampled animal as experimental unit. The Firmicutes to Bacteroidetes ratio was calculated for each group and analyzed with one-way ANOVA and orthogonal contrasts (i.e., CON vs MUL + MUH and MUL vs MUH). KEGG level 3 matrix was filtered for "KEGG level 2": this subset was then normalized and analyzed as described above. The significance level was set at 0.05. With respect to plasma and caecal metabolomics data, molecule concentrations were normalized via Box and Cox 76 transformation. Then, a one-way ANOVA and Tukey's post-hoc test were conducted considering the group as the experimental factor and each sampled animal as the experimental unit. To get an overview of metabolome trends, robust principal component analysis (rPCA) 77 was performed on molecules showing significantly different between-group concentration in the abovementioned univariate analysis. Initially, PcaHubert algorithm-implemented in rrcov package of R 67 -detected outlying samples according to their distance from others along and orthogonally to the PCA plane. Later, the optimal number of PCs was determined. A score plot and a correlation plot summarized main features of the rPCA models. The former represents the samples in the PC space, thus evidencing the overall structure of data. The latter reports Pearson's correlations between the concentration of each molecule and model components, thereby showing which molecule mostly affected the data structure. These analyses were carried out in R 67 with a significance level of 0.05.

Data availability
All data generated and analyzed during this study are included in this article and its additional files. Caecal microbiome sequences have been made public on MG-RAST repository with the project ID Muramidase_UniBO_ project_2020_34WGS (https:// www. mg-rast. org/ linkin. cgi? proje ct= mgp98 274).